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ABSTRACT 

We present one of the first physically-motivated two-dimensional general relativistic magnetohy- 
drodynamic (GRMHD) numerical simulations of a radiatively-cooled black-hole accretion disk. The 
fiducial simulation combines a total-energy-conserving formulation with a radiative cooling function, 
which includes bremsstrahlung, synchrotron, and Compton effects. By comparison with other sim- 
ulations we show that in optically thin advection-dominated accretion fiows, radiative cooling can 
significantly affect the structure, without necessarily leading to an optically thick, geometrically thin 
accretion disk. We further compare the results of our radiatively-cooled simulation to the predictions 
of a previously developed analytic model for such fiows. For the very low stress parameter and ac- 
cretion rate found in our simulated disk {a ~ 0.003, M/MEdd ^ 5 x 10~^), we closely match a state 
called the "transition" solution between an outer advection-dominated accretion fiow and what would 
be a magnetically-dominated accretion fiow (MDAF) in the interior. The qualitative and quantitative 
agreement between the numerical and analytic models is quite good, with only a few well-understood 
exceptions. According to the analytic model then, at significantly higher a or M, we would expect a 
fuh MDAF to form. 

The collection of simulations in this work also provide important data for interpreting other numer- 
ical results in the literature, as they span the most common treatments of thermodynamics, including 
simulations evolving: 1) the internal energy only; 2) the internal energy plus an explicit cooling func- 
tion; 3) the total energy without cooling; and 4) total energy including cooling. We find that the 
total energy formulation is a necessary prerequisite for proper treatment of radiative cooling in MRI 
accretion fiows, as the internal energy formulation produces a large unphysical numerical cooling of 
its own. We also find that the relativistic cooling functions must be handled carefully numerically in 
order to avoid equally unphysical heating or cooling runaways. 

Subject headings: accretion, accretion disks — black hole physics — galaxies: active — MHD — 
relativity — X-rays: stars 



1. INTRODUCTION 

The process by which turbulent accretion fiows around 
black holes develop large scale magnetic fields that can 
drive collimated jet outfiows is still poorly understood. 
Much of what we know comes from observations of X- 



ray binaries (XRBs) (summarized by Fender et al.|2 004), 
largely because the rapid variability in these stellar- 
mass systems allow for observations of state changes 
on timescales of days to at most a few years, whereas 
state changes are rarely observed for AGN. In XRBs, 
jets are asso ciated with both the Hard and Soft accre- 
tion states (McClintock & Remillard 2006), although 
their characteristics in the two states are quite differ- 
ent. Hard state jets are very steady, potentially last- 
ing for weeks, whereas jets produced in the transition 
to the high accretion-rate Soft state can be explosive 
and short-liv ed. One possible int erpretation of this phe- 
nomenology (Fender et al. 2004) is that accreting black 
holes produce jets most or the time, with the jet speed 
being a function of the accretion rate. Objects in the 
Hard state produce slow {v ~ 0.3c) jets, while Soft-state 
objects produce fast jets with fiat-space Lorentz factors 



W = (\ — v^l(?)~^l'^ ~ 10 that are comparable with jets 
in higher luminosity AGN (e.g., FR H radio sources). 
The explosive jet, in this picture, is simply the bow shock 
of a new fast jet interacting with a previously-existing 
slow jet as the source's accretion rate temporarily, and 
rapidly, increases toward the Eddington limit. 

Recent numerical simulations of non-radiative magne- 
tohydrodynamic (MHD) fiows that are unstable to the 
magneto-rotational instabi lity (MRI) have been shown 
to produce jetted outfiows ( [Hirose et al.|2004{ McKinney 
2006). In those simulations the mechanism involves the 
development of a magnetically-dominated region close to 
the black hole and rotation axis, where the magnetic field 
orders itself into a helical, rotating structure that drives 
the jet. At the present time, however, it is not clear if 
these results fit the phenomenology described above. The 
problem is that non-radiative MRI simulations should be 
the proper theoretical counterpart to the Hard state, but 
the simulations produce relativistic jets instead of the 
slow type jets associated with the Hard stat e. 

One of the present authors has suggested ( Meier|2005 



2008 [2009) that, at moderately low accretion rates 



where turbulent, advection-dominated accretion fiows 
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or ADAFs should occur), the inflow inside a radius 
Ri - 100 TG = 100 GM/c^ should develop a black hole 
magnetosphere str ucture similar to the force-free ones 
studied recentl y by [Tomimatsu fc TakahashT (2001); Uz- 
|densky| ( |2004[ |2005p . In this picture, closed field lines 
connecting the disk at Ri with the event horizon could 
funnel ionized plasma toward the black hole, creating a 
magnetically-dominated accretion flow or MDAF. Open 
field lines anchored near on the other hand, could 
drive a jet outflow with a speed set by the dynamical 
timescale near Ri [- [GM/RiY^'^ - 0.1c]. A further 
strength of this model is that a relatively large magneto- 
sphere might help explain why quasi-periodic oscillations 
in the Hard state are observed in the Hertz range, rather 
than kHz. 

Our suggested mechanism for MD AF/magnetosphere 
formation is radiative cooling in the previously-supposed, 
radiatively-inefficient ADAF. Cooling would lower the 
plasma pressure and decrease the disk vertical scale 
height, both of which lead to a dramatic increase in the 
dominance of magnetic stresses (ratio of magnetic to gas 
pressure greater than unity). Strongly-magnetized plas- 
mas are much more stable to the MRI, leading to a de- 
crease in turbulence and a more ordered magnetic field. 
This is pr ecisely the same process that occurred in the 
Hirose et al. (2004) and McKmneyj (|2006^ simulations, 
but now at a radius of ~ 100 re instead of a few. 

A critical assumption in this picture, though, is that 
the entire plasma (electrons and ions) cools. This is in 
contrast to current accretion theory, which asserts that 
whenever the flow enters a hot. Hard X-ray state, the 
transfer of thermal energy between ions and electrons is 
inefficient, leading to a two-temperature, optically thin 
ADAF with cool electrons and hot ions. However, some 
theoretical and numerical work has suggested t his may 
not be the case (e.g. Begelman & Chiueh IQSSj Sharma 



et al.||2007| ). There is also some observational evidence 
that efficient energy transfer from ions to electrons must 
occur even when a black hole is in a Hard X-ray state. 
Monitoring of many black-hole candidate sources shows 
that they can be found in the Soft state at bolomet- 
ric luminosities lower than in the maximum Hard state. 
This appears to be especially true in Hard states where 
a strong, steady jet is produced (e.g ., the plateau state) . 
Sources at the top right of Fig. 7 of Fende r et al. (2004) 
(the "FBG diagram") are quite hard and yet quite lu- 
minous. As a concrete example, Cygnus X-1 produces 
90% or more as much bolo metric luminosity in th e Hard 



state as in its Soft state (IMcConnell et al. 2002). This 
is hardly "radiatively inefficient" accretion by any stan- 
dard. The only truly inefficient state might be the Qui- 
escent state, of which the black-hole candidate A0620-00 
and the Galactic center black hole Sgr A* may be exam- 
ples. 

Our goal in this paper is not to resolve the contro- 
versy of whether or not Hard state objects radiate effi- 
ciently, but rather to investigate how they will behave 
if they do. We proceed by performing general relativis- 
tic MHD simulations of MRI-unstable black hole accre- 
tion flows, with two key differences from previous in- 
vestigations: we include in the energy equation a plau- 
sible high-temperature cooling function that is relevant 
for such flows, and we assume that electrons and ions are 



sufficiently thermally coupled that cooling of the former 
also cools the latter, keeping ~ Tg. For completeness, 
we actually compare four different classes of numerical 
models: one that evolves internal energy without includ- 
ing cooling, one that evolves internal energy and includes 
cooling, one that conserves total energy but does not in- 
clude cooling, and one that conserves total energy and 
includes cooling. These four classes of models span most 
of the simulations that have been carried out to date 
by other authors, and expands significantly on what has 
been done thus far with simulations involving physical 
cooling mechanisms. Our paper is unique in that it gives 
the first direct comparison of all four using a single nu- 
merical code and which, we believe, is the first to inves- 
tigate numerically the triggering of state transitions in 
cooled black hole accretion flows. 

2. NUMERICAL METHODS 

This work is carrie d out using the Cosm os++ astro- 



physical MHD code (Anninos et al. 12005'). Cosmos++ 



includes several schemes tor solving the MHD equations 
including a traditional artificial viscosity (AV) scheme 
and a new extended artificial viscosity (eAV) method. 
The AV scheme is based on an internal-energy-evolving 
(entropy-conserving) scheme, whereas the eAV scheme is 
a hybrid dual energy scheme that solves both the inter- 
nal and total energy equations. The eAV scheme has the 
obvious advantage that it conserves total energy; it is 
also potentially more accurate than other fully conserva- 
tive schemes in tracking the internal energy of the gas 
because of the dual treatment. Artificial viscosity based 
schemes, which both of these are, further have the ad- 
vantage that they are simpler to deal with when it comes 
to including extra physics such as the radiative cooling 
being added in this work. Furthermore, the combination 
of AV and eAV methods allows us to directly compare, 
within a single numerical code, the effects of including 
realistic heating and cooling processes in the evolution 
of MRI turbulent accretion disks. 

Cosmos++ has options to solve the MHD equations 
in either a Newtonian or general relativistic framework. 
Here the general relativistic form is used. In writing our 
equations we use the standard notation in which four- 
and three-dimensional tensor quantities are represented 
by Greek and Latin indices, respectively, and repeated 
indices imply summation. The equations of mass con- 
servation, momentum conservation, and magnetic induc- 
tion, common to both numerical methods used in this 
work, have the form 

dtD^d,{DV') = {) , (1) 
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(3) 
(4) 



where g^i, is the 4-metric, g is the 4- metric determinant, 
^ = V—gu^ is the relativistic boost factor, D = Wp 
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is the generalized fluid density, = /vP is the trans- 
port velocity, = g^^u^ is the fluid 4- velocity, 5*^ = 
W{ph -\- 2Pb)u^ is the covariant momentum density, P 
is the fluid pressure, Q is the artificial viscosity used for 
shock capturing, V (without subscripts or superscripts) 
is the adiabatic index, and Cp are coefficients to deter- 
mine the strength of the hyperbolic and parabolic pieces 
of the divergence cleanser, and A(p, T, B) is the cool- 
ing function of a gas with density p, temperature T, tem- 
perature scale height and magnetic field strength 5, 
as described in detail in the next section. (With indices, 
r indicates the geometric connection coefficients of the 
metric.) There are two representations of the magnetic 
field in our equations: is the 4- vector of the magnetic 
field, which can be defined in terms of the dual of the 
Faraday tensor (5^ = ii^*F^^), and B' = W{B'-B^V') 
is the boosted magnetic field 3- vector, where B'^ is recov- 
ered from the orthogonality condition B^u^ = 0, 

B° = -J [go^B' + Qi^B^V) . (5) 

The magnetic pressure is Pb = ||5|p/87r = 
gjj^iyB^B^ /Stt. We have assumed an equation of state 
of the form P = {T — l)pe wi th e being the internal en- 
ergy. We use the scalar Q from Anninos et al.| ( |2QQ5| ) with 



2.0 and ki = 0.3. We fix the divergence cleanser 
coefficients to be Ch — c^^l^x^^^j ISi and = c/j,, where 
Ccfl = 0.5 is the Courant coefficient, Axmin is the min- 
imum covariant zone length, and At is the evolution 
timestep. 

Both computational schemes also solve the internal en- 
ergy equation in the form 

dtE + d.^EV) = -PdtW -{P^Q) d^iWV) (6) 
-^WA{p,T,H,B) , 

where E = We = Wpe is the generalized internal energy 
density. The temperature T in the cooling function is 
recovered from the internal energy of the gas using the 
ideal gas law 

T = (r - l)e/{kn) (7) 

where n = p/{pmH) is the number density of the gas 
and we use p = 1.69. 

Additionally the hybrid dual energy scheme solves the 
following total energy equation 

dt£ + di{£V') = T,''-di{F') (8) 



where the total energy £ is defined as 



2Pb) 




the curvature source term is 

s° = -^T'^0r% , (10) 

and the divergence flux contribution is defined as 
F' = V^ {{g'' - g''V^) ((P + Pb)S] + Q]) (11) 
1 



Ideally this total energy equation would be sufficient by 
itself. However, total energy schemes can run into trou- 
ble when recovering local values for the internal energy. 
Defining Sd as the non-thermal or "dynamical" compo- 
nent of the conserved energy, 

, 2PbW^ 
= ^=H ^^ + V-^ 
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we write 



E 



-gw 



(12) 



(13) 
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{B'B^ - B^B^V' 



{r-i)g''{V^r ' 

for the internal energy extracted from the conserved en- 
ergy field. The trouble arises when numerical truncation 
errors accumulate to the point that the sum of different 
physical contributions exceed the total energy {£d > £)- 
This can occur in kinematic or magnetic field dominated 
flows and in the vicinity of strong shocks. The problem 
of negative energy can be avoided rather simply by forc- 
ing a minimum threshold on E to guarantee positivity. 
However, such a floor value is clearly not an accurate rep- 
resentation of the internal energy. Here we can benefit 
from having evolved the internal energy independently. 
We choose to only use the internal energy extracted from 
the total energy field whenever E > 10~^S. This avoids 
corrupting the internal energy value with numerical trun- 
cation error. The low density, background gas can also 
create accuracy problems for the total energy scheme be- 
cause the density is occasionally reset to a numerical floor 
value. Therefore, we further require £ > 10~^fmax as a 
condition for replacing E with E. This effectively ex- 
cludes the background gas. Finally, to recover as much 
disk heating as possible we always use the larger of E or 
E, provided the above two conditions are met. 

We find that the cooling timestep, Atcooi = Ccfie/A, is 
generally much shorter than the MHD timestep required 
for stability in the fluid evolution, A^mho = CcflAx/F, 
where Ax and V are characteristic zone lengths and ve- 
locities, respectively. Therefore, to save computational 
resources, we subcycle the cooling calculation, updat- 
ing the energy E^ as well as the temperature T and 
scale height in each zone using only the cooling 
"source" term [final terms in equation ([7|] and a timestep 
Atcooi until a full MHD timestep is reached, i.e. until 

Xliir^''(^^cooi)i = AtMHD- Then we update the total 
energy £ and momentum Sj according to the final terms 
in equations ([9| and (f3|, respectively. After that we pro- 
ceed with the next MHD update for all other evolution 
terms using the normal timestep AtMHD- Occasionally 
we have to deal with cooling timesteps that are unrea- 
sonably small due to very low temperatures or very high 
cooling efficiencies, regimes well outside the normal lim- 
its. To prevent the code from getting hung up at these 
points, we restrict N^teps to be < 100. This limit is usu- 
ally applied only in regions of very low density or very 
low energy where proper treatment of the fluid is inher- 
ently difficult. 

3. COOLING FUNCTION 

Three cooling processes are treated in this work: 
bremsstrahlung, synchrotron, and the inverse- Compt on 
enhancement of each of these two. Gener ally, we im- 
plement the equations of Esin et al. (1996), with some 
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changes. Below we describe first the equations that 
we use when the radiation is opticahy thin and some 
modifications that are necessary in order for the coohng 
computations to work weh in our numerical simulations. 
We then describe the modifications necessary when the 
plasma becomes optically thick to these radiative pro- 
cesses. In the extremely optically thick limit, the treat- 
ment is essentially the diffusion approximation. 

3.1. Optically Thin Limit 



The total cooling rate for the optically thin gas is ( Esin 
eFaI][T996l ) 

q~ = ^br,C + ^s,C Qs : (14) 

where q^^ and q~ are the bremsstrahlung and syn- 
chrotron cooling terms, respectively, and 7^br,c and 7^s,c 
are Compton enhancement factors. The details of how to 
compu te the Compton enhancements are given in (Esm] 



et al. 



( [1996 ). Basically, 7?(i/) is a modified exponential 

function of the Compton parameter 

y = 4(T + 4T2) (re. +t2J 

where T = kTe/rrieC^ is the dimensionless electron tem- 
perature, and Tes is the electron scattering optical depth. 
iri{v) is limited to a maximum value of SkT/hv^ where h 
is Planck's constant. 7^br,c is found implicitly by inte- 
grating ri{v) dq^^/ dv over appropriate frequencies, and 
7^s,c is approximated as ri{vc)i where Vc is the criti- 
cal frequency below which the synchrotron emission be- 
comes self-absorbed. Note: while we implement both en- 
hancements, because synchrotron emission is dominant 
at temperatures where Comptonization becomes impor- 
tant, only ?7s,c is important in our simulations. 

Th e un-Comptonize d bremsstrahlung cooling rate 
fromlEsin et al. (|1996|) is 



^br 



q± 



where 



np{ne + n+) 

1.50 X 10- 
2.12 X 10" 

Ce = (^e+^+) X 

J 2.56 X 10- 
\3.42 X 10- 

q^=nen^ x 

'3.43 X 10- 
6.84 X 10- 



X 

22 (1 + 1.781T1-34) 
22T[ln(1.123T + 0.48) - 



22Ti-^(l + l.lT- 
22T[ln(1.123T) ^ 

22 _^ -j^ 7T^) 

22T[ln(1.123T)^ 



_ ^2 _ 

1.28] 



(15) 
(16) 

T < 1 
-1.5]T > 1 

(17) 

1.25T2-^)T < 1 
T > 1 



(18) 



T < 1 
1.24] T > 1 



in units of erg cm ^ s ^. These represent cooling due 



to electron- ion (17), positron- ion (17), electron-electron 
(18), positron-positron (18), and electron-positron (19) 
processes. Here Up = rie — n+ is the number density 
of protons and n+ is the number density of positron- 
electron pairs. One can determine the ratio U-^/rip = 
(n+/ne)/(l — n+/ne) needed to calculate some of these 
terms using the following expression 



Up 



2 X 10 



ln(1.12T + 1.3)_ 

-4X^/2 exp(-2/T)(l- 
4( 



(19) 

0.015T) T < 1 



(112/277r)a2(lnT)3(l + 0.058/T)-i T > 1 



where af is the fine structure constant. 

The un-Comptonized synchrotron rate is a sum of op- 
tically thick and thin emission 



27rkT 



Jo Juc 



es{-u)dv 



(20) 



where H is the temperature scale height. The critical 
frequency can be found by equating the optically thin 
and thick volume emissivities at Vn 



(21) 



and solving the above expression numerically. For an 
isotropic, full Maxwellian distribution of electrons and 
positrons, the optically thin volume emissivity is (Ma- 
hadevan et aL||1996[ |Esin et al.|[T996l ) 



e,(z/,^)=4.43 X 10-''M7ri/(ne + n+) x (22) 
/(xM/sini?) _3 _i 

where xm = ^/^m is the normalized frequency (with 
vm = 6.27 X 10^^ B {kT^ [cgs] being the critical electron 
frequency for a given temperature), ^ is the angle be- 
tween the observer and the magnetic field direction, and 
K2 is the modified Bessel function of the second kind of 
order 2, given by the integral 



'Y^2 roo 

^ Ji/r 



' dz 



(23) 



In the high-temperature limit, the electron-energy- 
integrated, unitless spectrum , is given by the well-known 
expression ( |Pacholczyk|p7Q| 



/ ( ^) ^ ^ ^2 ^-z ^) (24) 

\smi9J Xm Jo 

where F{x) is the normalized synchrotron spectrum for 
a single electron 

F{X) = X / X5/3(0^^e 

J X 



|Esin et al. ( 1996| ) further average equation (23) over d to 
obtain the total emissivity needed in equations (l20|) and 
I2TI 



e^(iy)=4.43 X 10"^''47ri/(ne + n+) x 

I'{xm) _3 -1 

ergs cm s , 



(25) 



i^2(l/T) 



The angle- and energy-integrated, unitles s spectrum 
I'{xm ) can be fit to the following expression (iMahadevan 



et al.|1996; 



, 4.0505 / 0.40 0.5316 \ 

^(^m) = ^7^ 1 + ^ + ^7^ X (26) 



1/6 



1/4 



1/2 



exp(-1.8899x]^^) 

with no more than 2.7% error over the range < xm < 
00. 
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3.2. Problems with the Cooling Functions 



Generally, equations ( 14 



21) and (26 



K 



27) work fairly 
which is the 



well in the temperature range 10 
range over which they were used by |Esin et al.| (Il996| . 
However, this is not sufficiently broad for our numeri- 
cal simulations where the plasma temporarily can attain 
very high or very low temperatures in a given cell. In 
applying these equations over the temperature range ex- 
perienced in our simulations, we discovered the following 
problems: 

1. The number of positrons and electrons (determined 
from equation 20) diverges for T > 2.4 x 10^^ K, 
causing the sinmlation to crash. 

2. There is an error in the synchrotron cooling expres- 
sion (equation 26 ) that causes unphysical enhance- 
ment of the emission below T < 10^ K. The error 
is so severe, that without a fix the simulations de- 
velop a cooling runaway, which freezes the plasma 
into a cold, toroidally-dominated magnetic state. 

Our current fix for the first problem is very simple: 
we ignore positron cooling entirely (i.e., n+ = 0). A 
complete fix to the positron/electron ratio calculation is 
being investigated at present. However, our cooled sim- 
ulation generally remains below 10^^ K in most places, 
so our neglecting the contribution of positron cooling is 
reasonably valid. 

The second problem requires a more sophisticated so- 
lution. The error in the synchrotron cooling is caused 
by the use of different lower integration limits in equa- 
tions (23) and (24). i^2(l/'^) is the correct factor only 
if lower temperatures (T ^ 0) are allowed in equation 
(24) also. The result of this limit mismatch is that the 
denominator in equation ([26| vanishes faster for low tem- 
peratures (T < 10^ K) than the numerator, leading to 
enormous cooling rates for plausible temperatures (in the 
range 10^~^ K). S uch a problem would not have aflFected 
|Esin et ah] (1996)'s results, which maintained tempera- 
tures above this range. However, in a numerical simu- 
lation with many millions of cooling computations over 
millions of cells and time steps, the probability is quite 
high that such cool temperatures will be attained some- 
where in the flow, whereupon the entire structure will 
catastrophically freeze. 



Since we still wish to use equation (27) for I'{xm)^ the 
best fix for this problem is simply to assume the same 
high temperature limit in equation (23) as was assumed 
in equation ([24| (i.e., allow 1/T ^ 0, which is equivalent 
to replacing ir2(l/T) 2T^), resulting in a corrected 
synchrotron emissivity expression 

6^ (i/) =4.43 X 10"^° 2iTu {ue + n+) x (27) 

I'{xm) _3 -1 



to be used instead of equation rt26k in equation (20). 



The error introduced by using the high-temperature limit 
of K2 is of the same order as that caused by using a 
zero lower limit in the numerator of equation ( [24| ) (i.e., 
I{xm/ sin'i^)) in the first place. And that error is negli- 
gible compared to the total plasma emission, because it 
occurs at low temperatures where bremsstrahlung dom- 
inates (see Fig. [T]). 



Finally, a Saha ionization equation is used to deter- 
mine the electron density, which provides an exponential 
cutoff in the cooling below T ~ 10^ K. Thus, we consider 
continuum cooling only; no line or molecular cooling is 
included. 

One cautionary comment about these cooling functions 
is worth noting. The expressions contain many expo- 
nentials whose arguments easily can trigger underflow or 
overflow in a digital computer. If great care is not taken 
in respecting these limitations, even correct coding of the 
cooling functions will lead to non-physical results (heat- 
ing or cooling runaways) in the accretion flow. 

3.3. Optically Thick Limit 

To represent the cooling behavior in the opti cally 
thick l imit, we use a slightly modified version of |Esin| 
(1996)'s equation (21) that is suitable for multi- 



et al. 



dimensional MHD simulations. Our total cooling func- 
tion is given by Huben}^ ( 1990| ) 



A: 



-q 



2 ^^abs 



1 + V^Tabs 

AaT'^IH 

|r + \/3 + l/Tabs 



(28) 



(29) 



where the local temperature scale height is computed 
from 

This is a suitable definition for scale height in a multi- 
dimensional numerical simulation. If the flow were to 



assume a thin disk structure, for example, equation (30) 



would give the standard exponential (radiation energy 
density) scale height. 
The optical depth due to absorption is calculated as 



with 



'^abs — 



^abs 



(31) 



(32) 



And the total optical depth is computed using an aver 



aged opacity 

T = (k) pH 
with the diffusion average of k being 

i4M2 



-2 |V(T' 



pT^V • [V(T4)/(«p)] 
pH^V ■ [V{T^)/{kp)] 



(33) 

(34) 
(35) 



The total opacity used in this equation is given by 
^ = ^abs + i^es with /^es = 0.4 cm^ g~^ being the electron 
scattering opacity. Note that our definitions for H and 
(tv) allow for redistribution of heat within an optically 
thick region (r >> 1); equation (29) then reduces ex- 
actly to the diffusion approximation. Our approach also 
allows for non-local heating outside a marginally thick- 
thin transition region (r ^ 1) by photons from more op- 
tically thick reg ions, since the partially diffusive nature 
of equation ( 35 ) in this situation allows partial transport 



of heat from warmer to cooler regions. However, in the 
optically thin case (r << 1 and A = g~), there is no 
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Fig. 1. — Plot of a sample cooling function with p = 10 
g cm-3, B = 8380 G, and if = 2.7 x lO'^ cm. The total cool- 
ing function (solid line) and the following components are rep- 
resented: bremsstrahlung (short- dashed); Compton enhancement 
to bremsstrahlung (dot-long dash); synchrotron (dotted); and the 
Compton enhancement to synchrotron (dot- dot- dash). 

transport of heat from one region of the simulation to 
the other. There is only heat loss from the plasma and, 
therefore, from the simulation grid. 

For the purpose of illustration, in Figure [l] we plot the 
total cooling function A as a function of temperature, as- 
suming fixed values of p, and B. Of particular note 
is the temperature range over which each of the cooling 
processes is important. For 10^ < 10^ K, the domi- 
nant process is bremsstrahlung. At higher temperatures 
bremsstrahlung no longer dominates, even though its de- 
pendence with temperature steepens to be proportional 
to TlogT. Instead, above 10^ K the dominant cooling 
process is synchrotron, with Compton enhancement of 
synchrotron becoming important for temperatures a lit- 
tle above that. Comptonization of bremsstrahlung, while 
included, is never particularly important in our simula- 
tions. 

4. INITIALIZATION 

We initialize these simulations starting from the ana- 
lytic solution for a constant specific an gular momentum 
torus around a non-rotating black hole (Kozlowski et al. 
1978). In our initialization, the torus is defined by its 
mner radius rin = ISOrc and the radius of the pressure 
maximum reenter = 200rG. Knowledge of reenter leads 
directly to a determination of ^, the specific angular mo- 
mentum of the torus, by setting it equal to the geodesic 
value at that radius. Having chosen rin we can obtain 
Uin = Ut{rin)^ the surface binding energy of the torus, 
from Ut'^ = g^^ + i^g"^ 



The solution of the torus variables can now be 



speci- 



fied. The internal energy of the torus is (Hawley et al. 
T984| 

1 



(36) 



Thus the initial temperature of the torus, Tq = (F — 
l){limH /k)e^ is fixed to be ^ 10^ K by the specification 
of the torus. Assuming an isentropic equation of state 
P = pe(r — 1) = i^p^ for the initialization, the density is 



= 5 X 10^^ (cgs units). This gives an initial density 
maximum in the torus of Pmax,o = 2.8 x 10~^ g cm~^. 
Finally, the angular velocity of the fiuid is specified by 



(37) 



Once the torus is constructed, it is seeded with a weak 
dipole magnetic field in the form of poloidal loops along 
the isobaric contours within the torus. The initial mag- 
netic field vector potential is ( |De Villiers fc Hawley|2QQ3] ) 

A , - i Pent) for p > Pent , /ggx 

for p < Pent . ^ ^ 







The non-zero spatial magnetic field components are then 
— —deA(f) and = drA^p. The parameter pcut = 
0.5 * Pmax,o is used to keep the field a suitable distance 
inside the surface of the torus. Using the constant b in 
equation (38), the field is normalized such that initially 



P = P/P^ > Pq = 10 throughout the torus. The choice 
of the initial field geometry has been shown to have rel- 
atively little effect on the development of the M RI and 



the evolution of the disk (iBeckwith et al. 20081), which 

Tic 



given by p = [e{T - l)/f^] 



i/(r-i) 



We take F = 5/3 and 



is all we are focused on in this manuscript. However, 
the initial field topology does imprint itself in the forma- 
tion and evolution of jets, meaning that we will need to 
perform a more widely varying set of simulations before 
addressing that topic. 

In the background region not specified by the torus so- 
lution, we set up a static, low density (p = 10~^Pmax,o)5 
non-magnetic, hot gas. Numerical fioors are placed on 
p and e at approximately 10~^^ and 10 of their ini- 
tial maxima, respectively. The density fioor is very sel- 
dom applied once the initial background is replaced by 
evolved disk material. The energy fioor is applied some- 
what more frequently. Nevertheless, these very low fioor 
values should not have any significant dynamical impact 
on the problem. 

These simulations are performed in 2.5 spatial dimen- 
sions (all three spatial components of vector quantities 
are evolved, although symmetry is assumed in the az- 
imuthal direction) using a spherical polar coordinate 
grid. The grid used in the majority of the simulations 
consists of 192 radial zones and 128 zones in 0. We also 
performed select simulations at one-half and at double 
this resolution to test the numerical convergence of our 
results. We find very little variation between our default 
resolution and the higher resolution simulation, suggest- 
ing our results are well converged. 

In the radial direction we use a logarithmic coordinate 
of the form = 1.0 + ln(r/rBH)- The spatial resolution 
near the black hole horizon is Ar ~ 0.05rG; near the 
initial pressure maximum of the torus, the resolution is 
Ar ~ 5rG. Both are considerably smaller than the initial 
characteristic MRI wavelength Amri = 27rvA/^ ~ 50rG. 
In the angular direction, we use a concentrated latitude 
coordinate X2 of the form = X2 ^ ^{l — /i)sin(2x2) 
with h = 0.5, which concentrates resolution toward the 
midplane of the disk. As a result reenter A6> = Arc near 
the midplane while it is a factor of ~ 3 larger for the 
zones near the pole. 

For this work we have run the Cosmos++ numerical 
code in four different modes: 1) internal-energy evolv- 
ing with no explicit cooling (model 5221 or simply I); 
2) internal-energy evolving including an explicit cooling 
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TABLE 1 
Models 



Name 


Description 


I 


Internal-energy evolving 


I+C 


internal energy + cooling 


T 


Total-energy conserving 


T+C 


Total energy + cooling 




TABLE 2 




Parameters 


Name 


Initial value 



Tin 


ISOrG = 2.2 X 10^ cm 


'^center 


2mrG = 3.0 X 10^ cm 


^orb 


1.77 X lO^M = 0.875 s 


Pmax,0 


2.8 X 10-9 g cm-3 


^background 


2.8 X 10-1^ g cm-3 


^disk 


^ 10^ K 


^background 


^ 10^1 K 


Mbh 


1OM0 


MEdd 


8.4 X 10^8 g s-i 



function (model 522IC or simply I+C); 3) total-energy 
conserving with no explicit cooling (model 522T or sim- 
ply T); and 4) total-energy conserving including an ex- 
plicit cooling function (models 522TC or simply T+C). 
The motivation for this is to allow for a clear, direct com- 
parison of simulations carried out under different physi- 
cal assumptions. The "522" in the long naming conven- 
tion is a reference to our choice of = 5 x 10^^. 

5. RESULTS 

Since no cooling processes are treated in simulations I 
and T, those results simply scale with the mass of the 
black hole. However, for purposes of comparison with 
simulations I+C and T+C, we will assume the same scale 
for all variables in each simulation. Specifically we as- 
sume a black hole mass of M = lOM©, which sets the fol- 
lowing physical scales in the initial torus: rin = 2.2 x 10^ 
cm and reenter = 3.0 x 10^ cm. The orbital period at 
r = reenter IS ^orb = 1-77 X lO'^M = 0.875 s. The initial 
gas densities and temperatures are Pmax,o = 2.8 x 10~^ 

g Cm-3, pbackground = 2.8 X 10"^^ g Cm"^, Tdisk ~ 10^ 

K, and Tbackground ^ 10^^ K. The mass accretion rate is 
scaled by the Eddington rate MEdd = 8-4 x 10^^ g s~^ 
for an M = lOM© black hole. The models and parame- 
ters are summarized in Tables [T] and [2j Each simulation 
is evolved for seven orbital periods. This is sufficient 
time for all four models to achieve approximate equilib- 
riums inside r ^ = 150rG. However, because these 
simulations are carried out in two dimensions, the anti- 
dynamo theorem prevents a true steady-state from being 
achieved, so these are only approximations of the true 
state. 

5.1. Internal- Energy Evolving with No Cooling 

In this simulation, which we designate 5221 or just I, 
we only evolve the internal energy equation (equation [t]), 
ignoring the cooling function (A = 0). This mode of evo- 
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Fig. 2. — Pseudo-color plots of log(T) with contours of log p. 
The upper-left panel is the final time dump of the internal energy- 
model I; the upper-right panel is the final time dump of the internal 
energy plus cooling model I+C; the lower-left panel is the final time 
dump of the total energy model T; and the lower-right panel is the 
final time dump of the total energy plus cooling model T+C. The 
density contours are at p = 0.005, 0.016, 0.05, 0.16, and 0.5pmax,o- 



lution has commonly been used in the pa st (e.g., De Vil- 
liers fc Hawley|20Q3 Anninos et al.|2005 ), particularly in 



codes derived from the pioneering work ofWilson|(IT972 
From a thermodynamics perspective, running the code 
in this mode is an interesting case study. Because total 
energy is not conserved, any kinetic or magnetic energy 
dissipated in the disk (except through shocks) is simply 
lost from the simulation. In a sense, though, this creates 
a sort of thermodynamic equilibrium, wherein cooling (in 
the sense of energy lost from the disk) exactly matches 
dissipative heating everywhere in the simulation. Thus, 
without explicitly treating heating or cooling, this sim- 
ulation actually mimics one that includes heating plus a 
highly efficient cooling process. Some caution is in or- 
der, though, in making such a statement. Some heating 
and cooling mechanisms are captured in equation ([7|, 
speciffcally shock heating (through the artiffcial viscos- 
ity term) and adiabatic heating and cooling. These may 
not be balanced in the same way they would in a simu- 
lation that rigorously treated both heating and cooling. 
Furthermore, this treatment implies rapid, efficient cool- 
ing throughout the computational domain, regardless of 
ph ysica l conditions. We will explore this point further 
in §5.4[ In the upper-left panel of Figure |2j we plot the 
final distribution of gas density and temperature for this 
model. 

5.2. Internal- Energy Evolving with Cooling 

For this simulation, which we designate 522IC or sim- 
ply I+C, we again evolve the internal-energy equation 
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(equation [7|), this time including the radiative coohng 
term. Physicahy speaking, there is relatively little mo- 
tivation for this model, as we know there are important 
dissipative heating processes in disks that are ignored 
in this model. Nevertheless, this model does serve to 
round out our small lattice of tests and demonstrate the 
importance of using a fully conservative energy scheme 
(or some other heat-capturing procedure) when includ- 
ing radiative cooling processes. This is because, without 
including heating, there is nothing to counterbalance the 
cooling, and the disk ends up unreasonably cold and thin, 
with temperatures below 10^ K over much of the disk 
midplane, as shown in the upper-right panel of Figure [2j 

5.3. Total-Energy Conservation with No Cooling 

For this simulation, which we designate 522T or sim- 
ply T, we use the total energy conserving mode of Cos- 
mos++ (again with A = 0). By evolving equation (|9| 
and conserving total energy, we effectively capture dissi- 
pative heating mechanisms ignored in the previous simu- 
lations since any losses to the kinetic or magnetic energy 
of the gas are recovered as heat. Total energy conserving 
codes (as us ed previously by e.g., l McKinney|2Q'Q6| [Noble 
et al. 2007), are particularly applicable when consider- 



mg radiatively inefficient accretion flows (RIAFs), such 
as the one that is thought to be currentl y feeding Sgr A* 
( [Narayan et al. 1995 Yuan et al. 2003| ). Because these 
disks (both simulated and real) are not able to radiate 
their heat away efficiently, they tend to be very hot and 
vertically thickened, as shown in the lower-left panel of 
Figure |2[ 

5.4. Total Energy Conservation with Cooling 

For this simulation, which we designate 522TC or 
T+C, we use the dual energy evolving mode of Cos- 
mos++ described in ^ By including equation (|9| and 
conserving total energy, we again effectively capture dis- 
sipative heating mechanisms in the disk. In addition to 
the total energy equation, we simultaneously evolve the 
internal energy equation ^ to ensure we recover rea- 
sonable values for the internal energy whenever the total 
energy budget is dominated by non-thermal components. 
This is particularly important for calculating the temper- 
ature of the gas, which is a crucial input into the cooling 
routine. This is the only simulation where the dissipative 
heating processes are balanced by a physically motivated 
local cooling function, as described in ^ As expected, 
this leads to an intermediate disk state between the hot, 
thickened RIAF state of simulation T and the unreal- 
istically cooled disk in simulation I+C. The results are 
shown in the lower-right panel of Figure [2| 

6. COMPARISON OF NUMERICAL MODELS 

Simply looking at Figure [2] and comparing the four 
models, we already note a number of qualitative differ- 
ences. Obviously the disk in model T, which is expected 
to capture heating appropriately but includes only adi- 
abatic cooling, is much hotter and thicker than any of 
the other three simulations. This is consistent with the 
expectations of a radiatively inefficient, two-temperature 
gas, where the ions are poorly coupled to the electrons. 
The opposite extreme is model I+C, which includes ra- 
diative cooling processes without capturing most of the 
real, physical heating in the disk. This leads to a very 



thin, cold disk solution, which could only apply in cases 
of very weak turbulence or very low ionization. 

More interesting is to compare models I, the internal- 
energy evolving model, and T+C, the total energy plus 
cooling model. As we said before, model I can be thought 
of as a radiatively efficient model, but an unphysical one 
where cooling equals heating practically everywhere in 
the flow. This gives a much cooler disk than in model T, 
but also one in which the temperature increases mono- 
tonically as gas moves radially inward through the disk 
(compressive heating becomes more important). This 
is in contrast to model T+C, which shows an approxi- 
mately constant or even slightly decreasing temperature 
for r < 150rG, due to the efficiency of Compton enhanced 
synchrotron radiation. Next we make a more quantita- 
tive comparison of the models. 

6.1. Angle-Averaged Properties of the Simulations 

First, we construct density- weighted spherical shell av- 
erages of the various disk properties. The formula we use 
is ^ 

{Q)A{r,t) = jlJ£QV^ded4>, (39) 

r27T 



where A = f^^ dO dc/) is the surface area of the 

shell. We also average these quantities over the flnal two 
orbital periods of the simulations, 5torb = ^min ^ t < 
^max = 7torb5 to negate any transient features. The time 
averages are deflned as 



{Q)t 



1 



Qdt 



(40) 



The numerical results for the internal-energy model I, 
the internal-energy plus cooling model I+C, the total- 
energy model T, and the total-energy plus cooling model 
T+C are shown in Figure [Sj We have also included the 
predictions for the transition state solution, which we 
discuss below (Section |7.1[ ). 

There are clear outliers among the various disk models. 
For instance, the total-energy conserving model T ex- 
hibits signiflcantly lower density, pressure, and azimuthal 
magnetic fleld in the inner regions than any of the other 
three simulations. This actually owes to its much lower 
accretion rate (shown in Figure |5|; material is just not 
moving through the disk very quicEly. This, coupled with 
the considerably larger thickness of model T, leads to 
very low density and pressure. The internal-energy plus 
cooling model I+C is an outlier in the other direction, be- 
ing an order of magnitude cooler and thinner than model 
T. Models I and T+C, on the other hand, representing 
unphysical and physical cooling, respectively, look very 
similar in many regards. In fact, the only notable excep- 
tions are in T, a, and p. We mentioned the difference 
in T above, which is owing to the efficiency of Compton- 
enhanced synchrotron cooling for T > 10^^ K, and will 
return to the differences in a and P below. 

6.2. Volume- Integrated Properties of the Simulations 

In Figure [6] we plot the total integrated internal and 
magnetic energies for each of our four classes of mod- 
els. In the magnetic energy we can see the characteristic 
growth of the magneto-rotational instability on an or- 
bital timescale, after which it saturates. After about 2-3 
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Fig. 3. — Main disk properties plotted as a function of radius 
for the internal-energy model I, the internal-energy plus cooling 
model I+C, the total-energy model T, and the total-energy plus 
cooling model T+C. The data have been time-averaged over the 
final two orbital periods of each simulation. P, T, a, and f3 are 
density-weighted averages. The thick solid line in each frame is the 
solution for the MDAF transition region from equations ( |43| ). 

orbital periods the magnetic energy begins to decay due 
to accretion into the black hole, advection off of the grid 
through the action of jets and winds, and also due to the 
Cowling anti-dynamo theorem (the magnetic field is not 
able to regenerate itself in two dimensions). 

Initially there is very little heating in the disk as the 
magneto-rotational instability has not had time to build 
up the turbulence and the fiow is mostly laminar. Thus, 
in models I (which fails to capture heating realistically), 
I+C (which neglects heating, yet includes cooling), and 
T-hC (which captures heating, yet also includes cooling), 
the disk initially begins to cool. Once the MRI really 
kicks in after about 1 orbit, models I, T, and T+C begin 
heating, although model I heats more slowly than the 
two total energy conserving models. Model I+C never 
shows significant heating, clearly demonstrating that the 
local cooling always dominates. Once heating begins, 
the non-radiative models I and T never stop heating, 
whereas cooling appears to catch up with heating in the 
radiatively-cooled total-energy model T+C after about 4 
orbits. 

The differences between the internal energy curves of 
models I and T in Figure [6] give some indication of the 
amount of energy simply lost from the simulation by 
model I, which uses only the internal energy formulation. 
This amount of energy is comparable to the amount of in- 
ternal energy initially contained in the simulation. Like- 
wise, the difference between the internal energy curves 
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Fig. 4. Fig. [3] continued. B^, B'^, , V , and Va are 
density- weighted averages. 
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Fig. 5. — Plot of sign(m) log(l + m/10~^) as a function of radius 
for the internal-energy model I, the internal-energy plus cooling 
model I+C, the total-energy model T, and the total-energy plus 
cooling model T+C, where m = M/MEdd- The data have been 
time-averaged over the final two orbital periods of each simulation. 
By this time all three simulations have achieved a reasonably steady 
inflow solution for r < ISOrG?. 

of models T and T+C say something about the level of 
cooling in the disk. The physical cooling in model T+C 
is of the same magnitude as the unphysical cooling of 
model I, but as we have already shown, the resulting 
disk structure has significant quantifiable differences. 

We note here that the internal energy in the lower- 
left panel of Figure Ig] (model T) continues to increase 
all the way to the end of the simulation. This suggests 
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Fig. 6. — Plots of total internal (circular symbols) and magnetic 
(square symbols) energies as functions of time for the internal- 
energy model I (top-left panel), the internal-energy plus cooling 
model I+C (top-right panel), the total-energy model T (bottom- 
left panel), and the total-energy plus cooling model T+C (bottom- 
right panel). The energy scales have been normalized to the initial 
internal energy in the simulations. 

that dissipative heating of the disk has not yet been fuhy 
quenched by anti-dynamo processes, so the coohng seen 
in model T+C is genuine. 

7. COMPARISON OF NUMERICAL AND ANALYTIC 
RESULTS 

Ultimately we would like to make direct comparisons 
between our numerical results and observations of black- 
hole accretion disks in the Hard state. In the meantime 
we can also compare our numerical results with appli- 
cable analytic work. An interesting comparison can be 
made between our radiatively-cooled model T+C and the 
MDAF model. The basic idea of the MDAF is that catas- 
trophic cooling in the inner region of the disk should 
cause the disk to collapse vertically and dramatically re- 
duce the thermal energy relative to the magnetic, such 



that p = Pgas/Ps becomes < 1 (Meier 



2005 



2008). It 



is clear from Figure |2] that the radiatively cooled disk 
(T+C) indeed is much thinner than the uncooled disk 
(T). However, it is apparent from Figure [7| where we 
plot /3 over much of the domain of the simulation that 
although f3 is significantly lower in some regions of the 
radiatively cooled disk relative to the uncooled disk (par- 
ticularly for r < lOOrc), its angle average is not less than 
unity at any radius. 

Apparently we did not achieve a fully magnetically 
dominated state. Nevertheless, we can make a quantita- 
tive comparison between our numerical model T+C and 
the predictions of the MDAF model. Prior to becoming 
magnetically dominated, the model predicts that the in- 
flow should pass through a "transitional" state, in which 
/3 decreases from its initially large value to of order unity. 
It is this transitional inflow solution, then, that we wish 
to compare with our numerical simulation. 

7.1. Analytic Theory of Transitional Flow 

Analytic development of MDAF theory begins with the 
development of a simple analytic model for the ADAF 
structure in the region where the ion and electron tem- 
peratures are definitely equal (i.e., r > 144 re). The 
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Fig. 7. — Pseudo-color plots of (3 with contours of log p. Panel 
a is the final time dump of the internal-energy model I; Panel b 
is the final time dump of the internal-energy plus cooling model 
I+C; Panel c is the final time dump of the total-energy model T; 
and Panel d is the final time dump of the total-energy plus cooling 
model T+C. The density contours are as in Fig. [2] 

simple AD AF model is constructed i n a manner simi- 
lar to the [Shakura fc Sunyaev ( 1973 ) a-model, except 
that cooling of the flow is performed by an advective 
term [Qadv ^ MP/(27rr^p)] instead of the usual radia- 
tive term. The Compton parameter in this optically thin 
flow 



y 



16X2' 



remains less than unity for r > Rq, where 



Ro = 2.75 X lO^a"^/^ 



mm^^^ cm. 



(41) 



is defined as the radius where y = 1 m. the ADAF, 
with m — M/Mq^ and m = M/MEdd- Compton cool- 
ing is unimportant until the inflow approaches this ra- 
dius. Note that the electron scattering optical depth 
^es = i^es P H uses the electron scattering opacity hZes and 
disk scale height H. Because we use spherical geometry, 
the scale height is equivalent to 

H = r sin O 

where O is the disk angular scale height, with most of 
the accretion flow occurring in a polar angle range of 
7r/2 - e < 6> < 7r/2 + 9. 

Inside Rq^ y must be > 1, Compton cooling becomes 
important, and the transitional flow begins. In fact, the 
analytic MDAF models assume that generic Compton 
cooling is dom inant, so y ^ 1 must be true for r < Rq 
( Shapiro et al.| jl976). This cooling decreases the plasma 
temperature m the transition region, and therefore the 
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disk scale height, which also enhances and rB^ by 
compression. Conservation of magnetic flux in the steady 
MHD inflow requires that 

B'^ocr-^H-^ 

This means that the magnetic viscosity stress parameter 
a no longer can remain uniform with radius 



Since H ex r^/^ in this solution, the magnetic stress 
relative to the pressure must increase as the inflow ap- 
proaches the black hole {a ex r~^/^). The end of tran- 
sitional flow, and beginning of true MDAF flow, begins 
inside the radius 



Ri 



2.75 X 10^ a^^^^mrfi^^^ cm. 



(42) 



where a{r) becomes unity. For model T+C, we find a ~ 
0.003, m = 10, and m 5 X 10"^, so Rq ^ lUrc and 
Ri ~ Svg- We therefore do not actually expect a full 
MDAF solution for this particular model. Instead we 
can make comparisons with the transition solution that 
applies between Ri and Rq. 

The analytic transitional flow model model predicts 
the following set of sca ling relations for the disk proper- 
ties (Meier 20051 |2Q08D : 



Pc-- 
Pc- 
Tc-- 
H- 

a{x) - 
p{x)-- 

rB"^-- 

yr __ 

Va-- 



:2.1 X 10-^ a-' m- 
-7 A X 10^^ a-^/^ m 
:2.65 X 10^ a^/^ m-^/^ K 



-3/2 g cm ^ 



m X 
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cm s 
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(43) 



where x = r/Grc at the disk midplane. 

Figures [3] and |4] above include the predictions for the 
MDAF transition region (equations 43) over the ap- 
propriate radial range, Ri < r < Rq. We find that 
the radiatively-cooled numerical simulation T+C fits the 
MDAF transition solution remarkably well, except in cer- 
tain circumstances. 

First of all. Figure [3] shows that the density, pressure, 
and a{r) parameter are fit not only qualitatively but also 
quantitatively by the analytic transitional flow model, at 
least out to Rq ~ 150 tg- The temperature structure, 
however, is a factor of 3 cooler than the analytic model 
and is not as constant with radius. This is due to the 
analytic model's assumption that generic Compton cool- 
ing dominates (i.e., ^ ~ 1) when, in fact, it is specifically 
Comptonized synchrotron cooling that is important in 
the simulations. The latter can have a slightly differ- 
ent value and depends additionally on the magnetic field 
strength, resulting in a slightly lower temperature that 



may not be constant with radius. This discrepancy also 
affects the disk scale height, which scales as ~ VT. 

In Figure [4j again, most properties are fit well by the 
analytic model except for two: B^ and rB^. In the nu- 
merical simulation, while the general magnitude of the 
magnetic field (seen in Va/c) fits fairly well (as does 
even the axial component B^)^ the distribution of the 
rest of the magnetic field into r and (j) components ap- 
pears reversed from the analytic predictions. That is, the 
predicted flux conservation does not take place. There 
are different possible reasons for this: 

• Any radial shear that could create B^ from rB^ 
is suppressed by the 2-D, axisymmetric nature of 
the simulations. 3-D simulations may show the ex- 
pected distribution of B^ and B^ in the transition 
region. 

• The natural state of magnetized accretion flow, 
even with cooling, may be like that of the RIAF 
models (e.g., T and I): always dominated by 
toroidal magnetic field. In this case, the predicted 
MDAF would not arise even with catastrophic cool- 
ing. 

Therefore, it will be important to compare such 2-D 
simulations with similar 3-D ones to see how the inner 
transitional and predicted MDAF flows evolve when a 
third free dimension is added. Indeed, both toroidally- 
dominated and radially-dominated flows may be possible 
in nature in this region, with a state transition between 
the two occurring from time-to-time. 

8. CONCLUSIONS 

Using the Cosmos++ code, we have performed two- 
dimensional general relativistic MHD simulations of 
MRI-unstable accretion flows around black holes, with 
the potential of bremsstrahlung, synchrotron, and Comp- 
ton cooling of the high-temperature inflow. In the pro- 
cess of implementing the radiative cooling processes in 
our code, we made the following observations: 



The cooling function in Esin et al. 
valid in the range 10^ K < T < 10^ 



( |1996D , while 
K, needs spe- 
cial attention and care in order to be valid outside 
that range and not lead to heating or freezing run- 
aways in the simulations. 

• If radiative cooling is to be added to MRI simula- 
tions, then the energy equation also must properly 
handle the "viscous" heating caused by reconnec- 
tion and dissipation inherent in the MRI turbu- 
lence. In the present era of moderate-resolution 
MRI simulations (where dissipation is caused by 
numerical effects) , this can be handled in one of two 
ways: perform total-energy-conserving simulations 
and compute the internal energy by subtracting the 
kinetic and magnetic energies from the total; or 
use an artificial resistivity term to resolve current 
sheets and allow energy lost through numerical re- 
connection to be recaptured as heat. Although the 
artificial resistivity technique has been used with 
good success in many Newtonian applications (e.g. 
INitta et al.||2QQH [Stone & Pringle 2001; iFragile] 
et al. 2005), it has onl y recently been te sted in a 
relativistic MHD code ( |KomissarQv||2007 ). 
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• Indeed, evolving internal, rather than total, en- 
ergy without an additional procedure for recapturing 
lost heat produces an unphysical numerical cooling 
which can rival, or exceed, true radiative coohng. 
Furthermore, even when the magnitude of coohng 
is comparable, the resulting disk structure is quite 
different for the internal-energy-only model. 

If our assumption of Tg ~ Ti is valid, then we obtain 
the following results pertaining to the astrophysics of ra- 
diatively cooled, magnetized accretion flows: 

• Model T+C confirms the "transitional flow" solu- 
tion, which is proposed to connect an outer ADAF- 
like flow with an inner magnetically- dominated 
flow, as a viable MHD accretion inflow state. 

• The accretion rate and magnetic viscosity param- 
eter {M/MEdd - 5 X 10"^, a ^ 0.003) that result 
from our choices of inputs are in a range that pro- 
duces a large transitional flow region, without lead- 
ing to a completely magnetically dominated state 
(i.e. a and l3~^ never exceed unity before the flow 
enters inside the last stable orbit). 

• Comparison of the numerically-computed transi- 
tional flow and our prior analytic models of this re- 
gion show remarkable qualitative and quantitative 
agreement. Exceptions are limited to the temper- 
ature structure of the analytic model (which used 
a cooling model much simpler than the numerical 
functions herein) and the radial vs. azimuthal mag- 
netic structure (which likely was affected by the 
limitations of 2-dimensional axisymmetric MHD). 

Further investigations into the development of a true 
MDAF solution will require additional 2-dimensional 
simulations to investigate inflows with smaller transi- 
tional and larger predicted true MDAF regions (i.e., with 
greater a and M), and new 3-dimensional simulations 
to study the effects of cooling on the ratio of the ra- 
dial to toroidal magnetic field components. The answers 
to these questions will determine whether or not radia- 
tive cooling ultimately can trigger the formation of large 
black hole coronae (MDAFs). These simulations may 
also help confirm that MDAFs can form moderate-speed 
jets as proposed in our introduction. 

One final point should be noted. These are some of 
the first MRI simulations in which radiative cooling has 
an important dynamical effect on the accretion inflow. 
It was difficult, therefore, to predict what the resulting 
accretion rate and a parameter would be. The resulting 
values here {a = 0.003, M = 5 x 10"^), did not turn 
out to be appropriate for a source in the upper right- 
hand portion of the FBG diagram, as we had originally 
intended. In fact, they probably are more appropriate for 
a source in the lower right-hand portion with a rather low 
accretion rate. The applicability of these simulations to 
such a source, or really any source, will depend primarily 
on the validity of the assumption that Tg = T^, as made 
in our cooling model. 
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